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Numerous applications all the way from biology and physics to economics depend on the density 
of first crossings over a boundary. Motivated by the lack of analytical tools for computing first- 
passage time densities (FPTDs) for complex problems, we propose a new simple method based on 
the Independent Interval Approximation (IIA). We generalise previous formulations of the IIA to 
handle non-smooth processes, and derive a closed form expression for the FPTD in Laplace and 2 - 
transform space for arbitrary boundary and starting points in one dimension. We focus on Markov 
processes for which the IIA is exact. To apply our equations, we calculate the FPTD in two cases: 
the Ornstein-Uhlenbeck process and the discrete time Brownian walk. Our results are in good 
agreement with Langevin dynamics simulations. 


I. INTRODUCTION 

When the electric potential between the interior and exterior of a neurone exceed a certain threshold, the neurone 
fires. After firing, the interior potential is abruptly reset to its rest value and the process starts over. How often 
it starts over depends on external stimuli (e.g. light and touch) and firing frequencies of neighbouring neurones. 
To better understand neurone firing, and ultimately how neurones work, researchers mm use stochastic models to 
calculate how long it takes for the interior potential to pass the firing threshold for the first time. 

Neurone dynamics is by far not the only case where first-passage problems arise. Such problems frequently occur 
in physics, chemistry, biology, ecology and economics 011 ] and is one of the reasons why first-passage problems are 
so heavily studied. But despite enormous interest there are surprisingly few cases where we know the probability 
distribution of first-passage times analytically. Most cases are for Markov processes. 

First-passage time densities for Markov processes mainly comes from two approaches. In the first approach, the 
so-called method of images, one solves the Fokker-Planck equation with absorbing boundaries 00. Even though 
conceptually simple, it is limited to symmetric problems such as when the absorbing boundary is at the bottom of 
a symmetric potential well. The second approach is renewal theory mm- It works for non-symmetric problems 
but lead often to expressions in Laplace-space that cannot be inverted analytically. Even though useful, both these 
approaches are in practice limited to simple problems. In fact, neither of them can provide the first-passage time 
density for a Brownian particle in a harmonic potential for a general boundary and starting point. Thus, in order to 
address more complex first-passage problems we need better analytical methods. 

Another class of useful methods have been developed to solve persistence problems. In persistence problems 
one wishes to know the probability S(t) that a stochastic variable remains below or above a boundary from the 
start up to some time point t. The first-passage time density p{t) is simply related to the persistence according to 
p{t) = —dS{t)/dt. To calculate the persistence, some researchers [8UTT] use methods that enumerate all trajectories 
with an even number of boundary crossings and calculate the probability for each trajectory. But apart from a 
few special cases, these crossing probabilities cannot be calculated exactly and approximations are needed. One 
approximation scheme that gained popularity is the Independent Interval Approximation (IIA) |1214I4] , which assumes 
that the length of time intervals between consecutive boundary crossings are independent. However, in its present 
formulation the HA assumes that the processes has a well defined velocity which means that it cannot deal with non¬ 
smooth processes, such as Brownian motion or discrete processes. To apply HA to those processes these shortcoming 
must be remedied. 

In this paper we generalise the HA to non-smooth processes and discrete time series. Starting with the discrete 
case we find a simple expression for the probability density of first-passage times to a boundary from a general 
starting point in ^-transform space. We then generalise our equations to encompass the continuous case and obtain 
a similar expression but now in Laplace transform space. The expression is based on return probability densities 
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FIG. 1: Discrete-time stochastic process x{n) (e.g. the position of a particle) as a function of the number of time steps n 
{t = nAt where At is time increment). We denote the time spent above the boundary B hj T 2 ,Ti,... and below the boundary 
by TijTa,.... (Top) Original process x{n). (Bottom) Approximate x{n) where at each crossing event we draw a new position 
from the overshoot distributions A±(A). Note that A+(A) and A_(A) may be different where ’+’ (’ —’) means that the process 
is above (below) the boundary. 


to the boundary and the probability that the stochastic variable is above the boundary at some time. To show 
the applicability of our results we study two examples: the Ornstein-Uhlenbeck process (i.e. Brownian motion in a 
harmonic potential) and the discrete time Brownian walk. But our method is much more general and can in principle 
be used for any Markovian process and as an approximation for non-Markovian dynamics. 


II. METHODS 

In this section we outline the IIA framework and derive an expression for the first-passage time density (FPTD) for 
continuous and discrete processes in one dimension. We denote the FPTD by pB{t \ xo) where t is time, Xq = x(t = 0) 
is the starting point of the process and x = B is the location of the absorbing boundary. In discrete time we let n be 
the number of time steps and t = nAt where At is the time increment. To better understand the IIA approach we 
develop the mathematics for discrete processes and then show how it is generalised to the continuous case. 

The IIA equations herein relates three core quantities. The psin \ ccq), the probability ujy(n) that x{n) is above B 
at the nth time step given that a;o < B, and the return probability density that x{n) returns to B after a B-crossing 
either from above, or from below, ■!/’-(^); after n steps. The quantities w>(n) and '4’±(n) are inputs to our 

framework which one needs to calculate on case by case basis. The probability w>(n) is in general simple to calculate. 
We find it by integrating the probability density P{x,n\xQ) of x{n)-. 

pOO 

w>(n) = / P(a;, n|a;o)(ia;. (1) 

Jb 

The return probability densities on the other hand are more complicated and needs to be discussed further. 

To better understand consider a discrete process that pass through B repeatedly (see Fig. top). The 

number of steps that x(n) remains below B is denoted by Ti, T 3 ,..., and above B hy T 2 ,T 4 ,.... The number of steps 
to the first arrival, Ti, is special because it depends on xq. The density of Ti is simply pb{Ti\xo). After the first 
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B-crossing at Ti, x(Ti) ends up at some distance Ai > 0 above B, rarely precisely on B (i.e. Ai = 0 ). To calculate 
the distribution of T2, we must consider the trajectory from B + Ai back across B. We denote the distribution of 
T2 by '0+(72 j Ai) where we assume that the length of T2 is independent on Ti. This is the core assumption of the 
IIA and is true for Markov processes. At T2, the process crossed B from above and is below B by A2. To find 
the number of steps until the next crossing, T3, we must consider the trajectory from B — A2 back across B. The 
distribution of T3 is '0-(73) A2). Repeating this pattern we find '0_|_(T2i, A2i_i) {i = 1 , 2 , 3 ,...) for trajectories above 
B and ' 0 _ (T2i+i, A2i) below B, where the A’s are random numbers drawn from the overshoot distributions A±(A) 
(See Fig. bottom). If A is small with respect to B — xq, overshooting the boundary by A will not significantly 
change our final results for psin \ xq). We may therefore average ^/’±(n, A) with respect to A±(A): 

nOO poo 

' 0 _('u) = / 'tp-{n, A)X-{A)dA and = / ' 0 +(n, A)A+(A)dA ( 2 ) 

Jo Jo 

where the B dependence enter through ip±{n,A). In Appendix we show simulation results for the overshoot 
distribution for the discrete Brownian walk which is given by 

A(A) = y/|ertc (A) . (3) 

Working with the averaged return probability densities 'ip±{n) instead of ' 4 >±{n, A) implies that we ignore fluctuations 
in A and approximate the original process x{n) by a clipped process. The dynamics of the clipped process is: when 
x{n) crosses B, draw A from A±(A), make a jump to R ± A, and continue (see Fig. bottom). The clipped process 
is obviously different from the true x{n) but simpler to handle analytically. But the difference is small. We show in 
Appendix [b| for the discrete Brownian walk that w>(n) for the clipped process is practically indistinguishable from 
the true one. Below we formulate the IIA equations based on the clipped process. 

Based on the clipped process we may calculate a;>(A) (t^r = NAt) in terms of the number of R-crossings. Our 
derivation below is the discrete time version of the derivation for the continuous time case in [El [13]. Note, however, 
that in ISKISj the quantity psit \ xq) does not appear, as thermal equilibrium is assumed initially. Let Pk{N) be the 
probability for a trajectory starting in xq < B and ends up above B ai In after k crossings. uj-y{N) is then the sum 
of all such trajectories with odd number of crossings 


^>(^)= E (4) 

To calculate pi{N)^ assume that the first up-cross occurred dii ni < N and that there is no down-cross between rii 
and N. Since n\ can be anywhere from 0 to N^ this gives 

N 

Pi{N) = PB{ni\xQ)Q+{N - ni) ( 5 ) 

ni=0 


where the probability of not crossing is 


Q±(n) = 1 - y] V’zbK). (6) 

n '—0 

To find psiN), assume that the first up-cross occurred at ni, the first down-cross occurred between ni and 712, that 
the second up-cross happened between n2 and n^, and no down-cross between and N. This gives 

N N—ni N—n 2 

P3{N) = PB{ni\xo) E '*/'+(’^2) E i^-{n3)Q+iN - ns). ( 7 ) 

ni =0 712—0 71^—0 

Continuing this for p 5 ( 7 V), pt{N), ... leads to 

N N— Til N— 7l2 ^~'^2k-2 

P2k-i{N) = Y PB{ni\xo) Y E E tp-{n2k-i)Q+{N - n2k-i)- 

Tli —0 712—0 713 —0 n2fc - l=0 


( 8 ) 
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Summing over all odd number of i3-crossings we get a;>(-/V) [Eq. @]. To solve Eq. Q for | xq), we take the 
z-transform (/(z) = carry out the resulting geometric series. After some algebra we obtain 

Pb{z\xo) = -a;>(z)g{z) where g{z) = -^-. (9) 

This equation relates the first-passage time density to the probability a;>(z) of being above the boundary, and to the 
return probability densities 4’+iz) and ip_(z). This constitutes one of our main results in this paper. 


A. Independent Interval Approximation in continuous time 

The proper limit of our clipped process to a continuous time process is the following. When B is reached from 
below, the trajectory makes a jump to S + e, where e is a small constant. As e —>■ 0 we approach the continuous case. 
The overshoot distributions A±(A) for this cases is a Dirac delta function, which leads to 

/ OO 

'0±(t, A)i5 (A — (B ± e)) dA. (10) 

-OO 

In Appendix]^ we show explicitly how the clipped Ornstein-Uhlenbeck process convergence to the continuous one as 
e gets smaller. 

To derive the IIA equations in the continuous case, we proceed in the same way as for the discrete case using pk (t ), 
but with sums in Eq. ([^ changed to integrals —>■ ^ J dt) as we let At —)• 0 and n —)■ oo, while maintaining 
t = nAt constant. If we take the Laplace transform (/(s) = dt) of the sum over pk{t) [Eq. we obtain 

a similar geometric series as before that leads to 

PBis\xo) = suj>{s)g{s) where g{s) = ^ 

l-V'+(s) 

III. RESULTS 

In this section we apply our main results, Eqs. 0 and 0, to two examples. The first example is the continuous 
Ornstein-Uhlenbeck process, and the second one is the discrete Brownian walk. We also show that our results lead 
to Kramers escape for a general Gaussian stationary process (Appendix [C|), and that our method is consistent with 
the method of images (Appendix . To test the validity of our theoretical results we compare them to Langevin 
dynamics simulations (see Appendixfor simulation details). 


A. Application 1: Ornstein-Uhlenbeck process 

The FPTD for the Ornstein-Uhlenbeck process (OUP) is inherently difficult to calculate explicitly (TS]. The one 
exception is the symmetric case when the absorbing boundary is at the bottom of the harmonic well which can be 
solved with the method of images [16]. However, this does not work for the general problem. Instead, several efforts 
focused on the renewal equation in Laplace space. But because the renewal equation cannot generally be inverted 
analytically [13 H] researchers used numerical inversion m and series expansion around poles [50] . The last example 
m currently holds the best analytical approximations for the FPTD in the field. But even though in principle exact, 
none of their expressions are on closed form and must be evaluated numerically. To work with those expressions one 
must specify at least one cut-off parameter (sometimes two) which in practise must be done by trail and error. We 
compare our results to [50] using their so-called integral representation (see appendix for explicit details). Our 
approach does not share the cut-off parameter problem because we start in another end. We rely on a particular 
functional form of ip^(t) that is asymptotically true for Markovian Gaussian Stationary Processes (GSPs); They decay 
exponentially. The return probability densities, together with the probability density P(x,t|a;o), that we know, yields 
our final expression. 
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1. Analytical predictions 

To calculate the return probability densities, imagine first a symmetric GSP where x(t) behaves in the same way 
below and above the boundary B = 0. For Markovian GSPs we know the probability that x{t) does not change sign 
during the time interval [0, t] (i.e. the persistence) given that the process started in the infinite past such that it is 
stationary at t = 0. It is given by [S] 

Qit) = — arcsin (e~^*) ~ 6“’'*. (12) 

TT ^ ' 

Because of symmetry, the return probability densities are in this case equal. Denoting them by ^(t) and using 

that '0(t) = —dQ{t)/dt (because Q{t) = 1 — J* tj;{t')dt', see Eq. (H), we obtain that tp{t) decays exponentially 

~ (13) 

To generalise this expression to the asymmetric case B 0 where ip-it) '0+(^)j follow [8l [21] and introduce 
crossing rates r± from above and below B m- This leads to 

Q±(t) ~ exp(-r±t) (14) 

and 

ip±{t) ^ r±exp {-r±t). (15) 

With V'lblO hand, we may calculate the FPTD for a general Markovian GSP from Eq. ([TT]). First, we use that 
i/'±(s) = \/{l + s/r±) in Laplace space. Second we put t/’±(s) in g[s) [Eq. 0] so that (/(s) = l + r+/(r_+s). After 
inversion we find 


PB{t I Xo) 


dlOy (t) 

dt 


+ r+ 


dt 


dt'. 


(16) 


This equation is valid for any Markovian GSP. To get further, we need to address specific examples where we can 
calculate r± and iOy,{t). 

For the OUP, a;> (t) is given by 


= Lrfc (^;£2LL) , 


It is calculated from P{x,t\xo)dx where |7] (dimensionless variables, see Appendix | a| 


P{x,t\xo) = 


•\/27r(l — e 2*) 


exp - 


{x-xpe *)^ 
2(1 - e-2‘) 


(17) 


(18) 


To calculate crossing rates r±, we proceed as follows. First, we get r_(_ from the normalisation condition 
f^ao Psit I ^o) dt = 1, or pbIs —>■ Ojxo) = 1. Also, using g{s —>■ 0) = 1 + r+jr- and lims_,,o sa;>(s) = w>(t —)■ oo) in 
Eq. ^ leads to 


r+ 


r. 


erfc [B / ^/2) 



(19) 


The other crossing rate r_ is equal to the inverse of the mean-first passage time t to B from xp = 0, where the 
latter result follows from the backward Fokker-Planck equation |22| 


r. 


1 

r 


where 


r = 




( 20 ) 


This can be understood as follows. When B is far away from the potential minimum x = 0, up-crossing events 
are rare. When they are rare, the distribution of times between up-crossing events is approximately the same as 
the FPTD. Because the process’s equilibration time is much shorter than the mean-first passage time to B, xp is 
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FIG. 2: First-passage time density ps(t | xq) out of a harmonic well V(x) = a;^/2 when the boundary is at B = 3 and xo = —1 
(left), Xo = 0 (middle), xo = +1 (right). The solid line is Eq. (161, ’o’ are results from Langevin dynamics simulations (averaged 
over 10® realisations), and ’x’ is the approximation in [20] [Eq. (Ellj. The insets show the behaviour at short times. 


approximately equilibrated with average xq = 0 [28] • This means that the FPDT is asymptotically Kramers expression 
PB(t|0) ~ [l^. This implies that 'ip(t) ~ psitjO) with r_ = 1/r. 

Before we compare Eq. (161 with simulations we clarify some of its limitations. First, note that Eq. is exact 
within the IIA but in general we do not know '!/’±(0- If assume that '<p±(t) decays exponentially we arrive at 
Eq. (16). This assumption is only asymptotically true. This means that Eq. ([T^ is less accurate when B and 


Xq are too close. If they are, several crossing events occur at short times and is not decaying as a single 

exponential. To better understand what is meant by ’short times’ we see in Fig. that B = 3 and xq < 0 gives 
good agreement with Langevin dynamics simulations for pB(t\xo), while there is a discrepancy for short times when 
B = 3 and xq > 1. To improve the results for short times we could in principle take into account sub-leading terms 
tp±{t) = r±e“’'=‘=‘[l + + ...]. However, this would introduce new parameters as well 


as conditions for when to truncate the sum. Using the asymptotic behaviour of 'fp±(t) makes our approach free of 
cut-off parameters. 


2. Simulations and numerical results 


To validate our method we compare pB(t|xo) to Langevin dynamics simulations. In Fig. [^we show pB{t\xo) for 
different Xq keeping B = 3 fixed. The simulation results are represented by circles and the IIA Eq. (16) by solid lines. 


The main panels show the overall behaviour whereas the insets show the short time dynamics. Overall, we get good 
agreement with simulations. 

In the long time limit, our analytical FPDT decays exponentially. This agrees with Kramers escape rate, and 
our IIA equations capture this regime well for all values of B — xq- For short times, there is a discrepancy between 
simulations and our IIA equations when B — xq gets smaller. This is because a considerable amount of first-passage 
events occur for short times before has attained its single exponential form, as discussed above. Also, the 

reason why our IIA equations systematically underestimates the simulations for short times is the following. Assume 
for simplicity thatxg = 0 and that B infinitesimally above Xg. In this case, the first-passage dynamics is practically 
indistinguishable from Brownian motion for short times. For Brownian motion we can show that ip±{t) is a Dirac 
delta function (Appendix which means that there is an infinite number of boundary crossings once the boundary 
is crossed. This is very different from ~ that says that the average time between two boundary crossings 

is l/r±. 

In Fig. [^we also included one of the best analytical approximations for the FPTD |5D] (see Appendix [^. Their 
formula approximates the short time dynamics better than our method while for long times both approaches match 
well with each other. 


3. First-passage through two boundaries 


It is straightforward to generalise our method to two boundaries, B > 0 and B' < 0. Clearly, ii B \B'\ we need 
four return probability densities which our method is unable to handle. But for the symmetric case, B = \B'\, ip±(t) 
are enough to describe the crossing in and out of the region —B<x<B. To calculate r_ we use a generalisation 
of Eq. (20) to two boundaries (see e.g. [22]), and to get r_|_ we use Eq. (19) where we replace a;>(oo) —>■ 2w>(oo) 
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FIG. 3: First-passage time density pB{t\xo) out of a harmonic well V{x) = a:^/2 with two boundaries and xo — 0: (left) 
B = ±2.5, (middle) B = ±3, (right) B = ±3.5. The solid line is Eq. (161, and ’o’ are results from Langevin dynamics 
simulations (averaged over 10® realisations). The insets show the behaviour at short times. 


(the probability that |a;(t)| > B). Figure]^ shows that our IIA equations match simulations for two boundaries with 
similar accuracy as for one boundary. 


B. Application 2: Discrete time Brownian walk 

The Brownian walk [3] is one of the simplest cases of a non-stationary process. Nevertheless, its FPTD is not 
known for general B and Xq except in terms of a double Laplace transform [9] that no one thus far have been able 
to invert. The one exact result that exists is the Sparre Andersen theorem [53] that says that the persistence to stay 
above (or below) the boundary when B = xq is 

( 21 ) 

Here we use q{n) and the IIA formalism to put forward a simple summation formula for the FPTD for general B and 
Xq. As before, we first need to find the return probability densities 'i/’j_(n). 


1. Analytical predictions 

In this problem x(n) behaves in the same way on both sides of the boundary, which means that the return probability 
densities on either side of B are equal, = tpin). We approximate 'ijj{n) with the discrete derivative of q{n)-, that 

is ipin) « — [q{n) — q{n — 1)] 0(n — 1). Here 0(n) is the unit step function (discrete heavy side step function) that 
takes care of the initial condition '0(n = 0) = 0. From Eq. (21) it follows that 

If we put Eq. © yields g{z) = 1 Using this in pBiz\xo) [Eq. ©] and inverting to n-space 

leads to 


PBin\xo) = w>(n) ± ^ w>(n - A:) [^(/c) - ip{k - 1) - 4,i] 


(23) 


k=l 


where = ^erfc [(B — xo)/(v^)] is calculated in the same way as Eq. 0- Equation ( |^ is a generalisation 

of the Sparre-Andersen theorem to general boundary and initial conditions. 

For long times we expect that PB{n\xQ) — V’(fi)- Indeed, expanding Eq. ( |22| ) for large n we get the Brownian walk 
result pB{n\xQ) ~ 


2. Simulations and numerical results 


In Fig. I^we compare Eq. (23) to simulations. Overall we find good correspondence, especially as B — xq increases. 
But as it decreases, we start to see deviations for small times, e.g for B — a:o = 1. The reason is that the overshooting 
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FIG. 4: First-passage time density pB{n\xQ = 0) for the simple Brownian walk when the boundary is at: B = 1 (left), 5 = 3 
(middle), B = 5 (right). Connected crosses come from Eq. (231 while connected rings represent simulations (averaged over 
2 X 10^ realisations). Insets display short-time dynamics. 


start to play a role and the derivative of the persistence is no longer a good approximation to ipiji). From simulations 
we find that the average overshooting length is 0.626... that is comparable to S — xq = 1. 


IV. SUMMARY AND OUTLOOK 

There are plenty of examples where one wants to know the probability density of first-passage times to a boundary. 
To find this density, one often use the method of images or renewal theory. These approaches are however in practise 
limited to simple cases. To find better methods, we improved the so-called Independent Interval Approximation (IIA), 
developed for survival probability problems, that is limited to smooth stochastic processes where trajectories have 
well defined velocities. This excludes for example Brownian motion. We generalised IIA to continuous and discrete 
non-smooth trajectories. From our IIA formalism we derive a simple expression for the first-passage time density 
to a general boundary and initial condition in one dimension. This expression relies on that we know the functional 
form of the return probability densities. But once it is known our approach is parameter free. To show the validity 
of our expression we applied it to the Ornstein-Uhlenbeck process (OUP) and the discrete time Brownian walk. For 
the OUP we use that the return probability densities decays exponentially for long times [5], ~ exp {—r±t), where we 
identify r_ (up-crossing rate) as the reciprocal of the mean first-passage time, which is known |22j . Then, from the 
normalisation condition of the first-passage time density we determine the last parameter, r_|_ (down-crossing rate). 
In discrete time we apply our IIA formula to the Brownian walk where we use the Sparre-Andersen theorem which 
yields the return probability densities. Both cases match well with Langevin dynamics simulations. We also show 
that (i) our approach reproduces Kramers expression for escape of a Brownian particle out of a harmonic potential, 
and (ii) that it becomes equivalent to the method of images for symmetric problems, e.g. when the boundary is at 
the bottom of a potential well. 

Our IIA equations are new, and we anticipate that they will have a wide applicability to previously intractable 
first-passage and escape problems. 
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Appendix A: Simulations 

Following [5S] we simulate the Ornstein-Uhlenbeck process with 

x(t -(- At) = x{t) -I- \fl 


_g-2Ai_^(0,l) 


(Al) 
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FIG. 5: (left) Convergence of the clipped process to the continuous Ornstein-Uhlenbeck process as e —>■ 0. For smaller and 
smaller e the simulated a;>(t) (Langevin dynamics, see Appendix]^ gets increasingly closer to the analytical result. Here we 
put B — 1 and xq = 0. (right) Simulated ojy{n) as a clipped process compared to the analytical result for the discrete Brownian 
walk with i? = 3 and xo = 0. (inset) Distribution of overshooting length A(A) when the boundary B — 0 is crossed for the 
discrete Brownian walk. All simulations are averaged over 10® realisations. 


where A/'(0,1) is a normally distributed number with mean zero and variance one. This equation is on dimensionless 
form where we made the replacements 


X —)■ 



and 


7 


(A2) 


where ksT is thermal energy, k the harmonic spring constant and 7 is hydrodynamic friction. The diffusion constant 
is D = ksT/^. In the simulations we varied xq and B and made statistics of when x(t) reached B for the first time. 
We averaged over 10® — 10^ ensembles. 

It is well known that the Langevin dynamics scheme systematically overestimates the first-passage time because it 
can potentially miss crossings that happened within At. To reduce this error we used adaptive time steps that get 
smaller as x(t) approaches B. We change At as follows: 

1. Set At = Ato. 

2. Calculate the probability that x(t + At) is above B given that x(t) is below B. That is w>(At). 

3. If a;>(At) > e then At —/Si 12. Otherwise do not change At. 

4. If at a later time ujy{At) < e, then At —>■ 2At with Ato as upper limit. 

In the simulations we used Ato = 10“® and e = 10“^. 


Appendix B: Supplementary figures 

In Fig. [^we illustrated how we approximate the real process x{n) by the so-called clipped process. Here in Fig. 
we show explicitly how the clipped process differs from the real x{n) in for the continuous Ornstein-Uhlenbeck process 
and the discrete time Brownian walk. To the left in Fig. [^we show the continuous time case. Clearly we approach 
the analytical curve for (t) [Eq. 0 ] as we let e —>■ 0 [Eq. 0 ] in our simulations. To the right in Fig. we show 
the discrete case. We found the overshoot distribution A(A) = ■\/ 7 r/ 2 erfc (A/2) from simulations (see inset). Using 
A(A), we simulate w>(n) as a clipped process. We see good agreement with analytics in all aspects. 
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x{t) 



Appendix C: Kramers escape 


For long times our method is consistent with Kramers escape theory. To see this we use the final value theorem 
and find that Eq. 0 in the limit s —>■ 0 (long times) becomes after Laplace inversion 


^ w>(oo)r+e 


Using Eq. (191 to eliminate r+ gives 


PB{t):^r_{l-ujy{oo))e 


(Cl) 

(C2) 


When B is large enough we get that 1 — w>(oo) « 1, Kramers expression [23ll26]. As an example, for the harmonic 
potential we find that w>(oo) « 0.001350 when B = 3 and Xq = 0. 


Appendix D: Method of Images 

The method of images can successfully give the EPTD to a boundary for symmetric Markovian problems. Here we 
show that our IIA formalisms is consistent with this method for Gaussian processes. To show this, we must first find 
the return probabilities. 

Eor a symmetric problem, e.g. when the absorbing boundary is a the bottom of a harmonic well, the return 
probability densities on either side of the boundary are equal, ip+{t) = ■*/'-(0 = '*/'(^)- To find them we first realise 
that the Langevin equation close to the bottom of the well, say a: = 0, is similar to unbiased Brownian motion: if 
the force is —x, the Langevin equation is dx{t)/dt = —x + ri{t) « ry(t), close to x = 0, where ry(t) is white noise. For 
Brownian motion we may use the following trick to find 

Consider a Brownian particle that diffuses between two boundaries separated by a small distance e (see Fig. |^. 
We denote an ’up-cross’ by crossing B from below and a ’down-cross’ by crossing B from above. When x = B is 
crossed the particle jumps immediately to i? ± e, as explained in Sec. |llj Now, the distribution of times between an 
up-cross and a down-cross is just the first-passage to a point that is a distance e away, that is [7] 

with Laplace transform 

?/’(s) = exp (-|e| yi) . (D2) 

If we now let e —>■ 0, up- and down-crossings occur to the same boundary and therefore = 1, or 

= (D3) 

This result manifests the fractal nature of Brownian motion: if there is one H-crossing at time t, there will be infinitely 
many in the infinitesimal interval (t, t + dt) |16j . 
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With ip{t) at hand we may derive the method of images formula. First we put tp^s) = 1 in Eq. ( |11[ ), which after 
after inversion leads to 


PB{t\xo) = 2 


dujy {t) 
dt 


(D4) 


Second, since pBit\xo) = —dS{t)/dt, S {t) i s the probability of not crossing a given boundary up to time t, and 
S'(O) = 1 and w>(0) = 0, we rewrite Eq. ( |D4[ ) as 

S{t) = l-2u}>(t). (D5) 

For a Gaussian process with mean and variance w>(t) is 


= ^erfc 


f B- p{t) \ 


(D6) 


Using Eq. 
formula. 


(D6| in (D5) and the relations erfc(*) = 1 — erf(*) and erf(—•) 


—erf(*) leads to the method of images 


S{t) 


1 

2 


1 + erf 




1 

2 


1 + erf 


/ p{t) -b \ 


[P{x, 1 1 xo) — P(x, 1 1 2B — ccq)] dx. 


(D7) 


Appendix E: Alili’s formula 


To compare our result for the Ornstein-Uhlenbeck process to one of the best known approximations we have 
implemented one of the formulas from |20j . In our notation it reads 


2t H_Aii2t){-BlV2) 

f H-AI(2t)-k-^ilt{-Xol V2) \ 
i , ^ \ A /+ ( — T3 / \/ 2 ^ / 


(El) 


where H^{z) is the Hermite function of order v. Here A and N are parameters that are determined based on trial and 
error. We found that for t < 10 —)■ {A = 18.1, N = 1000} while for t > 10 —)■ {A = 7,N = 1000}. The comparison to 
Langevin dynamics simulations and our IIA formula are seen in Fig. 
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